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We perform a principal component analysis (PCA) of vz ,{ pt ) in event-by-event hydrodynamic 
simulations of Pb-fPb collisions at the Large Hadron Collider (LHC). The PCA procedure identifies 
two dominant contributions to the two-particle correlation function, which together capture 99.9% 
of the squared variance. We find that the subleading flow (which is the largest source of flow 
factorization breaking in hydrodynamics) is predominantly a response to the radial excitations of a 
third-order eccentricity. We present a systematic study of the hydrodynamic response to these radial 
excitations in 2-l-lD viscous hydrodynamics. Finally, we construct a good geometrical predictor for 
the orientation angle and magnitude of the leading and subleading flows using two Fourier modes 
of the initial geometry. 


I. INTRODUCTION 


Two-particle correlation measurements in ultrarela- 
tivistic heavy ion collisions provide an extraordinarily de¬ 
tailed test of the hydrodynamic description of heavy ion 
events. Indeed, the measured two-particle correlations 
exhibit elliptic, triangular, and higher harmonic flows, 
which can be used to constrain the transport properties 
of the quark gluon plasma (QGP) produced in heavy ion 
collisions [1, 2]. In hydrodynamic simulations of heavy- 
ion events, fluctuations in the initial state are propagated 
by the expansion dynamics of the QGP, and this expan¬ 
sion ultimately induces fluctuations in the momentum 
spectra of the produced particles. Thus, measurements of 
the momentum space fluctuations (or correlations) con¬ 
strain the properties of the QGP expansion and the ini¬ 
tial state. The purpose of the current paper is to classify 
and quantify the dominant momentum space fluctuations 
in (boost-invariant) event-by-event hydrodynamics, and 
then to optimally correlate these fluctuations in momen¬ 
tum space with specific fluctuations in the initial state 
geometry. The current paper is focused on triangular 
flow, since it is a strong signal and driven entirely by 
fluctuations [3]. The corresponding studies of the other 
harmonics are postponed for future work. 


Due to flow fluctuations the correlation matrix of 
event-by-event triangular flows, (V3 (pti)V 3 *(pt 2)), in hy¬ 
drodynamics does not factorize [4]. Factorization break¬ 
ing is quantified by the parameter r{pTi,PT 2 )-, 


r{pTl,PT2) 


{VsjPTiWi {PT2)) ^ 

VmiPTlW) {\V3{PT2W) ~ ’ 


( 1 ) 


which must be less than unity when there are several 
statistically independent sources of triangular flow in the 
event sample [4]. Factorization breaking has been stud¬ 
ied in event-by-event hydrodynamics [4-6] and compares 


reasonably to the measured data for appropriate param¬ 
eters [6]. It is generally understood from these analy¬ 
ses that factorization breaking is caused by the hydro- 
dynamic response to geometrical properties of the initial 
state that are poorly characterized by the coarse geomet¬ 
rical measure £3. For instance, in Ref. [6] the r{pTi,PT 2 ) 
matrix was found to be sensitive to a parameter control¬ 
ling the roughness of the initial state. In Ref. [5] it was 
suggested that a careful study of the r matrix and other 
observables could be used to test hydrodynamic predic¬ 
tions for the pt dependence of the event plane angle, 
which arises when multiple triangular flows are present 
in a single event. The current paper clarifies the origin 
of factorization breaking by associating the largest non- 
factorizable contribution to the triangular flow with the 
hydrodynamic response to the first radial excitation in 
the triangular geometry. 

First, in Sec. II we use principal component analysis 
(PCA) of the harmonic spectrum to analyze the trans¬ 
verse momentum dependence of the third harmonic in 
boost invariant event-by-event hydrodynamics. PCA is 
a statistical technique that decomposes the flow corre¬ 
lation matrix into eigenvectors and eigenvalues [7]. The 
procedure naturally identifies the most important contri¬ 
butions to flow fluctuations. Typically only two modes 
are needed to give an excellent description of the full 
covariance matrix {V^ipriW^{ jpt 2 )) to 0.1% accuracy. 
When there are only two significant eigenvectors (or tri¬ 
angular flow patterns), the r matrix can be expressed 
as [7] 


r{pTi,PT2) - 1 “ 2 


/ vt\pTi) 

\V3^\pti) 


V3^\pt2) '\ 

V3^\pT2)) 


( 2 ) 


where V^^\pt) and V^‘^\pt) are the first and second 
eigenvectors^. 


^ As described in Sec. II, the eigenvectors are normalized to the 
eigenvalue and we are assuming that 

Ai ^ A2. 
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The leading mode of the third harmonic is strongly 
correlated with the triangular event plane[3] , and thus is 
essentially equivalent to familiar measurements of v^Ipt) 
with the scalar product or event plane method. How¬ 
ever, the subleading mode is uncorrelated with the lead¬ 
ing event plane, and is therefore projected out in most 
measurements of harmonic flow. Section III studies the 
basic properties of the subleading triangular flow, such 
as its dependence on centrality and viscosity. 

In Sec. IV A we show that the subleading triangular 
flow arises (predominantly) from the radial excitation of 
the triangular geometry. To reach this conclusion we 
first directly calculate the average geometry in the event 
plane of the leading and subleading flows. This averaged 
geometry (as explained in Sec. IV A) is shown in Fig. 4 
and exhibits a familiar triangular shape for the leading 
flow and a triangular shape with a radial excitation for 
the subleading flow. 

Having identihed the physical origin of the sublead¬ 
ing flow, we introduce several geometric predictors which 
(with various degrees of accuracy) quantitatively predict 
the magnitude and orientation of the subleading flow in 
event-by-event hydrodynamics based on the initial data, 
in much the way that £3 predicts the orientation and 
magnitude of the leading vz{pt)- 

As a first step, in Sec. IV B we correlate the principal 
momentum space fluctuations with the Fourier modes of 
the geometry. Based on this analysis in Sec. IV C we 
construct a good geometrical predictor for the orientation 
angle and magnitudes of the leading and subleading flows 
based on two Fourier modes. 

For comparison, we also correlate the subleading flow 
with a linear combination of the complex £3 3 and £ 3 , 5 , 


£ 3.3 - n3 

(3a) 

^3.5 - ^5 • 

(3b) 


'rms 


where the square brackets [] denote an average over the 
initial entropy density in a specific event, and i?rms = 
([r^]) is the event averaged root-mean-square radius. 
Note that our definitions of £ 3,3 and £ 3,5 are chosen to 
make the event-by-event quantities £ 3,3 and £ 3^5 linear 
in the fluctuations, since the denominator is a constant 
event-averaged quantity. In this respect this definition is 
different from the conventional one which is a nonlinear 
function of the initial perturbations^ (see Sec. HB for 
further explanation). We find that the subleading mode 
is also reasonably correlated with a linear combination of 
these two quantities, but the quality of this predictor is 


^ We compared analogous results with £3 3 and £3 5 defined via 
cumulants [8], and found them marginally worse than the ones 
presented in this paper. 


considerably worse than a predictor based on two specihc 
Fourier modes. 

The geometric predictors described above are ulti¬ 
mately based on the assumption of linear response. 
At least for the third harmonic (the scope of this pa¬ 
per), these assumptions are checked in Sec. V. In this 
section we explicitly compare the response to the av¬ 
erage (“single-shot” hydrodynamics [9]) and the aver¬ 
age response (event-by-event hydrodynamics). We find 
reasonable agreement between these two computational 
strategies for both the leading and subleading triangular 
modes. 


II. PCA OF TRIANGULAR FLOW IN 
EVENT-BY-EVENT HYDRODYNAMICS 


A. Principal components 


PCA was recently introduced in Ref. [7] (which in¬ 
cludes one of the authors) to quantify the dominant mo¬ 
mentum space fluctuations of harmonic flows in trans¬ 
verse momentum and rapidity in a precise way. This sec¬ 
tion provides a brief review of this statistical technique. 

Paraphrasing Ref. [7] , in the flow picture of heavy ion 
collisions the particles in each event are drawn indepen¬ 
dently from a single particle distribution which fluctuates 
from event to event. The event-by-event single particle 
distribution is expanded in a Fourier series 

rlN °° 

^ = ^o(PT) + E^-(P^)e”“^ + H-C-, (4) 

^ n—1 


where dp = dy dpr dp notates the phase space, p is the 
azimuthal angle of the distribution, and H.c. denotes 
Hermitian conjugate. Vnipr) is a complex Fourier coeffi¬ 
cient recording the magnitude and orientation of the nth 
harmonic flow. This definition deviates from the common 
practice of normalizing the complex Fourier coefficient by 
the multiplicity, n„(pr) = K(pt)/4o(pt)- 

Up to non-flow corrections of order the multiplicity N, 
the long-range part of the two-particle correlation func¬ 
tion is determined by the statistics of the event-by-event 
fluctuations of the single distribution 


dApairs 

dpidp2 


dN dN 
dpi dp 2 


+ 0{N) . 


( 5 ) 


If the two-particle correlation function is also expanded 
in a Fourier series 




-in{ipi-ip2) 


,dpidp2, 

then this series determines the statistics of 14 ,(pt) 


( 6 ) 


VnA{PTl,PT2) = {Vn{PTl)V*{pT2)) ■ (7) 
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The covariance matrix t 4 ,A(PTi,PT 2 ), which is real, sym¬ 
metric, and positive-semidefinite, can be decomposed 
into real eigenvectors, 

Kia(PT1,PT2) (8) 

a 

= (9) 

a 

where Vn°‘\pT) = (pt) and /“dpr = 

6 ab- As discussed above we have not normalized V 2 ,(pt) 
by the multiplicity. To make contact with previous work, 
we define and present numerical results for 

||»Wf ^^, (10) 

f {dN/dpr) dpT f {dN/dpr) dpx 


plane according to a two-component model. Specifically, 
for the ith participant we assign a weight 


Ai = K 


(1 — a) a 

2 ^ -^{n-coiiji 


(13) 


with a = 0.11, K = 35.1 for rj/s = 0.08, and k = 32.8 
for p/s = 0.16. (ncou)i is the number of binary col¬ 
lisions experienced by the ith participant; so the total 
number of binary collisions is A"coii = | X)i(ncoii)i- The 
entropy density in the transverse plane at initial time Tq 
and transverse position x = {x,y) is taken to be 


s{To,x)= ^ S^{To,X-Xi), (14) 

^^-^parts 


where Xi = (x, y) labels the transverse coordinates of the 
ith participant, and 


which scales with multiplicity and £3 in the same way as 
an integrated ^ 3 ( 2 } measurement. Typically in event-by¬ 
event hydrodynamics (as shown below) the eigenvalues 
are strongly ordered, and two eigenvectors describe the 
variance in the harmonic flow to 0.1% accuracy. Thus, 
PCA provides a remarkably economical description of the 
momentum dependence of flow fluctuations. 

The harmonic flow in each event can be decomposed 
into its principal directions, 

VsIpt) = ^iV3^^\pt) + ^ 2 Vi^\pT) + ■■■■ (11) 

The real vectors V^^\pt),V^‘^'^ (pr), ■ ■ ■ (which do not 
fluctuate from event to event) record the root-mean- 
square amplitude of the leading and subleading flows. 
The complex coefficients , ^ 2 , ■ • ■ indicate the orienta¬ 
tion and event-by-event amplitude of their respective 
flows. The amplitudes of the different components are 
uncorrelated by construction 

{Ut)=Sab. ( 12 ) 

The original impetus for this work was a desire to under¬ 
stand which aspects of the geometry are responsible for 
the orientation angle of the second principal component. 


B. Simulations 

In this paper we use boost-invariant event-by-event hy¬ 
drodynamics to study the principal components of V 3 (jjt) 
for LHC initial conditions. The implementation details 
of the hydrodynamics code will be reported elsewhere, 
and here we note only the most important features. Our 
simulations are boost invariant and implement second or¬ 
der viscous hydrodynamics [ 10 ], using a code base which 
has been developed previously [11, 12]. For the initial 
conditions we use the Phobos Glauber Monte Carlo [13], 
and we distribute the entropy density in the transverse 


1 _ 

s^{to,x) = A ,--—^ (15) 

with '/2<t = 0.7 fm. The parameters k and a are 
marginally different from Qiu’s thesis [14], and we have 
independently verified that this choice of parameters re¬ 
produces the average multiplicity in the event. ^ 

The equation of state is motivated by lattice QCD 
calculations [15] and has been used previously by Ro- 
matschke and Luzum [16]. In this paper we compute 
“direct” pions (i.e. pions calculated directly from the 
freeze-out surface) and we do not include resonance de¬ 
cays. We use a freeze-out temperature of Tfo = 140 MeV. 

Simulation results were generated for fourteen 5% cen¬ 
trality classes with impact parameter up to b = 12.4 fm 
and at two viscosities, p/s = 0.08 and p/s = 0.16. Un¬ 
less specified, the results are for p/s = 0.08. We gen¬ 
erated 5000 events per centrality class. ^ We then per¬ 
formed PCA for the third harmonic Vs^px) by discretiz¬ 
ing V 3 {pt) results from hydrodynamics into 100 equally 
spaced bins between px = 0 ... 5 GeV, and finding the 
eigenvalues and eigenvectors of the resulting Hermitian 
matrix. Similar results for the other harmonics will be 
discussed elsewhere. 

Table I records the Glauber data which is used in this 
analysis. Event-by-event averages with the initial en¬ 
tropy density are notated with square brackets, e.g. 



® More precisely we have verified that for these parame¬ 
ters hydrodynamics with averaged initial conditions reproduces 
^ function of centrality after all resonance de¬ 
cays are included. Assuming that the ratio of the charged parti¬ 
cle yield to the direct pion yield is the same as in the averaged 
simulations, the current event-by-event simulations reproduces 
dNch/dr). 

^ We thank Soumya Mohapatra for collaboration during the initial 
stages of this project. 
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Centrality 

(^min 5 

^max) 

^ part 

-^rms 

^ max 

_rms 

^3,3 

0-5 % 

(0.0, 

3.3) 

384 

4.1 

8.1 

0.11 

5-10 % 

(3.3, 

4.7) 

335 

3.9 

7.8 

0.14 

10-15% 

(4.7, 

5.7) 

290 

3.7 

7.5 

0.17 

15-20% 

(5.7, 

6.6) 

250 

3.6 

7.3 

0.20 

20-25% 

(6.6, 

7.4) 

215 

3.4 

7.0 

0.22 

25-30% 

(7.4, 

8.1) 

184 

3.3 

6.7 

0.25 

30-35% 

(8.1, 

8.8) 

156 

3.2 

6.4 

0.28 

35-40% 

(8.8, 

9.4) 

132 

3.1 

6.2 

0.32 

40-45% 

(9.4, 

9.9) 

no 

3.0 

5.9 

0.35 

45-50% 

(9.9, 

10.5) 

91 

2.9 

5.7 

0.39 

50-55% 

(10.5, 

11.0) 

74 

2.7 

5.4 

0.44 

55-60% 

(11.0, 

11.5) 

60 

2.7 

5.1 

0.48 

60-65% 

(11.5, 

11.9) 

47 

2.6 

4.8 

0.52 

65-70% 

(11.9, 

12.4) 

37 

2.4 

4.4 

0.55 


TABLE 1. Table of parameters from the Glauber model (all 
distances are measured in fm). 


where S'tot is the average total entropy in a given cen¬ 
trality class, (/ (Px ToS{to, a:)). Averages over events are 
notated with (), so that the root mean square radius is 

(17) 

As a technical note, here and below the radius is mea¬ 
sured from the center of entropy, so [a:] = 0 . £ 3^3 and £3 5 
are defined in a somewhat unorthodox fashion in Eq. (3), 
with £™ = a/(| e 3 , 3 p)- is the averaged maximum 

participant radius, max|a;i|. 


III. SUBLEADING TRIANGULAR FLOW 


As a first step, we list the (scaled) magnitudes of flows 
11 ^ 3 “^ II [Eq. (10)] in central collisions for the simulations 
described above: 


a 

1 

2 

3 

4 

^ CO 

1.5 X 10-2 

2.6 X 10-3 

4.8 X 10-4 

1.1 X 10-4 


Note that the quantities in this table are proportional to 
the square-root of the eigenvalues, 11 ^ 3 “^ || oc From 

the decreasing magnitudes of the listed (scaled) magni¬ 
tudes, we see that the first two eigenmodes account for 
99.9% of the squared variance, which can be represented 
as a sum of the eigenvalues 


we divide by (dN/dpT) and present the same eigenmodes 
in Fig. 1(b). 

We also investigated the centrality and viscosity de¬ 
pendence of the principal components. The normalized 
principal flow eigenvectors ^“(pr) are approximately in¬ 
dependent of viscosity (not shown). In Fig. 2, we show 
the centrality dependence of these normalized eigenvec¬ 
tors. In more central collisions the eigenvectors shift to 
larger transverse momentum, which can be understood 
with the system size scaling introduced in Ref. [17]. 

The magnitude of the flow, i.e. the squared integral 
/ {y^°'\pT)YdpT, depends on both centrality and vis¬ 
cosity. To factor out the trivial multiplicity dependence 
of VsIpt), we plot the scaled flow eigenvalues ||ui“^|| [see 
Eq. (10)] in Fig. 3. Going from rj/s = 0.08 to p/s = 0.16 
we see significant suppression of the leading mode. In 
general the subleading scaled flow depends weakly 

on centrality. 


IV. GEOMETRIG PREDICTORS FOR 
SUBLEADING FLOW 


A. Average geometry in the subleading plane 


In this section, we clarify the physical origin of the sub¬ 
leading flow by correlating the subleading hydrodynamic 
response with the geometry. 

As a first step, we determined the average initial geom¬ 
etry in the principal component plane. Specifically, for 
each event the phase of the principal component [see 
Eq. (11)] defines orientation of the flow. We then rotate 
each event into plane and average the initial entropy 
density, S{x) = Tos{to,x). More precisely, the event- 
by-event geometry in the principal component plane is 
defined to be 


Si.r, (j); ^a) = \ S (r, (j) + {aig^a + 2Tr£)/3) , (19) 

where we have averaged over the phases of Next, 

we average S{r,(p; ^a) over all events weighted by the 
magnitude of the flow 

5(r,<^;ea) = (5(r,<^;Ca)|ea|). (20) 


dpT {Vz{pt)V;{pt)) Aa oc ^ ||4“^f . (18) 

a a 

Figure 1(a) displays the eigenvectors, V^°'\pt), for the 
leading and first two subleading modes. We see that 
only the first two flow modes are significant, and in the 
rest of this paper we consider only these two. To make 
contact with the more traditional definitions of vs^px). 



Figure 4 shows the in-plane averaged geometry S{r, 4>] ^a) 
for the leading and subleading principal components in 
central collisions. Clearly, the leading principal compo¬ 
nent Ug is strongly correlated with the triangular com¬ 
ponents of the initial geometry, while the subleading com- 
ponent Eg is correlated with the radial excitations of 
this geometry. 

To give a one-dimensional projection of Fig. 4, we in- 
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FIG. 1. Momentum dependence of flow components in central collisions, a) Principal flow vectors, b) Principal flow 

vectors divided by the average multiplicity, v^‘^\pt) = V^°‘\pt)/ {dN/dpr). 



Pt (GeV) 

FIG. 2. Gentrality dependence of flow eigenvectors 'i/)“(pT). 
tegrate Eq. (20) over the azimuthal angle to define 
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subleading x 
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ijxxxxxxxxxxx 

(^5) IppHinn J 


FIG. 3. Gentrality and viscosity dependence of scaled eigen¬ 
values 11 ^ 3 “^ II. (The subleading flow has been magnified 5 
times to bring to scale with leading flow.) 


p2TT 

S3{r-,^a)= / ( 21 ) 

This is equivalent to defining 53 ( 7 -), 


, 53 ( 0 = / dct>S{r,cj,)e^<^, ( 22 ) 

and correlating this with the flow fluctuation 

Mr-, ia) = {SM)Ca) ■ (23) 

Results for S' 3 (r; Mr^ are shown by the blue (gray) 
curves in Fig. 5. 


Again we see that the leading flow originates from a 
geometric fluctuation with a large integrated eccentric¬ 
ity, while the subleading flow is sensitive to the radial 
excitation of the triangularity. Note that the relatively 
small subleading flow corresponds to a fairly significant 
fluctuation of the initial geometry. 
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FIG. 4. Average geometry xr^ in the leading and subleading principal component planes in central collisions minus an averaged 
radially symmetric background, r^{S{x\ ^a) — (S(a:)|^a|)). Peak fluctuations are ±10-20% above the background. 



FIG. 5. Gorrelation between the principal components and 
the triangular geometry, (S 3 (r)^*), for the leading and sub¬ 
leading flows in central collisions. The result has been mul¬ 
tiplied by and normalized by StotR^ms, so that the area 
under the leading curve is approximately £ 3 ™°. 


Here S^ik) has the meaning of the mth harmonic of the 
2D Fourier transform of the initial geometry S'(fe), i.e. 

We recall that the £ 3,3 is determined by the long wave¬ 
length limit of S 3 {k) [ 8 ] 

r o ri \ ~o {kRy-^s/2)^ 

hm Sslk) = -Stot --£ 3,3 , (25) 


where i?rms is the root mean square radius and S'tot is 
the total entropy in a given centrality bin. The constant 
factors are determined by the expansion of Jsix) near 
a; = 0. Motivated by this limit we define a generalized 
eccentricity e 3 (fc) 


(26) 

which approaches £ 3,3 as fc —f 0. Clearly in a Glauber 
model there is an analogous definition 


esik) = / 

‘^tot J 0 


rdrSs^r) 


Wo, 




30i 


part 


3! 


(fci?,.^,/2)3 


Jsikvi) 


, (27) 


B. The average geometry in Fourier space 


It is evident from Fig. 5 that the leading and sublead¬ 
ing geometries have different characteristic wave num¬ 
bers. This becomes apparent when we correlate the flow 
signal with the Fourier (or Hankel) transform of the tri¬ 
angular geometry 


Ssik) ^ f 


rdr Sz{T)J^{kr) . 


(24) 


where the coordinates of the fth participant are Xi = 
(r cos tpi, r sin (pi). 

The Pearson correlation coefficient between the flow 
and a specific wave number k is 


Qa{k) 


{L4ik)) 

x/(k3(fc)P)(|5aP)' 


(28) 


Examining Qa{k) in Fig. 6 , we see that leading compo¬ 
nent is produced by low-/c fluctuations, while subleading 
flow originates from fluctuations at larger k. 
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FIG. 6. Quality plot (or Pearson correlation coefficient) for 
e^ik) as a single-term predictor for principal flows (central 
collisions). 


C. Optimal geometric predictors for the subleading 

flow 


In this section, our aim is to predict the magnitude 
and orientation of the leading and subleading flows. To 
this end we regress each principal component of the flow 
with various Fourier components of the initial geometry. 

Following Ref. [18], we construct a prediction for the 
flow amplitude by taking a linear combination of 

£z{k) : 


^pred _ E uJbSsih) ■ (29) 

i=l 

The selected wave numbers kb are discussed in the next 
paragraph. The response coefficients tOb are chosen to 
minimize the square error , or equivalently to maximize 
the Pearson correlation coefficient Qa between the flow 
and the prediction 

min ( 30 ) 

/c c*pred\ 

max Qa = /- . ( 31 ) 

The correlation coefficient Qa is referred to as the quality 
coefficient in Ref. [18]. 

We construct two predictors based on two and five 
wave numbers. For the five-term predictor we choose 
equidistant points which span the range seen in Fig. 6 

fcbiirms = 1, 3, 5, 7, 9, (32) 

and fit the response coefficients uJb- The two term pre¬ 
dictor was motivated by the discrete Fourier-Bessel series 


advocated for in Ref. [19], 

kbRo J3,l5 J3,2 77(3 — 3 Ri'nis 5 (33) 


where (j 3 ,i,j 3 , 2 ) — (6.38,9.76) are the zeros of Jaix), 
and we select Ro to optimize the correlation between the 
geometrical predictor and the flow. For comparison, we 
also constructed a two term linear predictor from the 
familiar eccentricities £ 3^3 and £ 3^5 defined in Eqs. (3a) 
and (3b). 

In Figs. 7(a) and (b) we show the correlation coefficient 
between the leading and subleading flow amplitudes ^a 
and the predicted amplitudes using the two and five 
Fourier mode fits, and the £ 3 ^ 3 , £ 3^5 fit. As is well known, 
the leading mode is very well predicted by £ 3,3 and £ 3 ^ 5 , 
though the quality degrades towards peripheral collisions. 
These results for the leading mode can be compared prof¬ 
itably with Fig. 4 of Ref. [20], where similar results were 
recently reported. For the subleading flow the linear cor¬ 
relation coefficient is reduced relative to the leading flow, 
and a high degree of correlation is only achieved for the 
0-40% centrality range. The simple geometric predictor 
based on £ 3^3 and £ 3^5 is reasonably correlated with the 
subleading flow in central collisions, but this correlation 
rapidly deteriorates in more peripheral collisions. The 
predictor based on two judiciously chosen wave numbers 
generally outperforms all other predictors we studied for 
both the leading and subleading modes. Indeed, we be¬ 
lieve based on numerous other fits that these two wave 
numbers essentially exhaust the predictive power of the 
event-by-event triangular geometry, S^^r). Additional 
histograms correlating the amplitudes and phases of the 
flow and the two-term prediction are shown in Fig. 12 in 
the appendix. 

The two-wave-number fit correlates the flow with a 
specific projection of the triangular geometry, i.e. 

pOO 

rdrSsiQpir), (34) 

Jo 

where p{r) is a radial weight chosen to maximize the 
correlation between the flow and the projection. This is 
analogous to using £ 3^3 to predict triangular flow, where 
the radial weight is p{r) oc — 


pOO 

£ 3.3 0 ^ - / rdrS3{r)r^ . 
Jo 


We have used Fourier modes as a basis for p{r) 


p{r) oc wi 


Joikir) 


0 J 2 - 


Jo{k2r) 


{kiR 

rms )3 (fcai? rms 


(35) 


(36) 


but other functions could have been used.^ In Fig. 8 we 
compare the radial weights for the leading and subleading 


^ A table of ujiIuj2 is given as a function of centrality in the ap¬ 
pendix. 
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FIG. 7. Pearson correlation coefficient for 2 and 5 fc-mode predictors [Eq. (26)], and a predictor based on the eccentricities £ 3,3 
and £ 3 , 5 . (a) Correlations for the leading flow (zero suppressed for clarity), (b) Correlations for the subleading flow. 
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FIG. 8 . Radial weight functions for the leading and sublead¬ 
ing flow predictors in central collisions, see Eq. (34). 
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ometry, we formed the analogous correlator between the 
predicted flow and 5'3(r), (-S' 3 (r)^*P"'®'^). Examining 
Fig. 9, we see that the two term predictor fully captures 
the optimal average geometry. For peripheral collisions 
the optimal geometry differs from what we can construct 
using linear combinations of Fourier modes, suggesting 
that additional nonlinear physics [9, 12, 18] plays a role 
in determining the subleading flow. 

Figure 7(b) also shows the correlation (or lack thereof) 
between the subleading flow and the integrated V 3 


(6^3*) 


(38) 


Since the subleading mode is uncorrelated with the lead¬ 
ing mode (by construction), there is almost no correla¬ 
tion between the integrated V 3 and the subleading mode. 
The upshot is that measurements of V 3 {pt) based on the 
event plane or scalar product method are projecting out 
the important physics of the subleading mode. 


modes. The overall normalization of weight function p(r) 
is adjusted so that 

rdrp{r)S^{r) 

The weight function for the leading projector is very close 
to cubic weight, but the subleading radial weight has a 
node at r ~ 1.5i?i.ms • 

Within the framework of linear response, in Sec. IV A 
we found the optimal geometry for predicting the lead¬ 
ing and subleading flows by correlating the observed flow 
with the geometry, (S' 3 (r)^*). To test if the two and 
five wave number predictors reproduce this optimal ge- 



- 


(37) 


V. TESTING LINEAR RESPONSE 


The success of the linear flow predictors discussed in 
previous section depends on the applicability of linear re¬ 
sponse. A straightforward way to check this assumption 
is to compare the averaged response of event-by-event hy¬ 
drodynamics to the hydrodynamic response to suitably 
averaged initial conditions. 

In Sec. IV A we computed the average geometry in 
the event planes of the leading and subleading flows (see 
Fig. 4). It is straightforward to simulate this smooth 
initial condition and to compute the associated V 3 (pt)- 
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FIG. 9. Comparison of the averaged geometry in the principal 
component plane, {S 3 {r)^a), and two- and five-term predictor 
plane, in central collisions. 



0 0.5 1 1.5 2 2.5 3 3.5 

PT (GeV) 


FIG. 10. Comparison of event-by-event hydro (averaged re¬ 
sponse) and single-shot hydrodynamics (response to average 
geometry) in central collisions. The singe-shot hydrodynamic 
results are generated from the initial conditions in Fig. 4. 


This is known as “single-shot” hydrodynamics in the lit¬ 
erature [9]. In Fig. 10 we compare V^{pt) from the lead¬ 
ing and subleading average geometries to the principal 
components V^^\pt) and V^‘^\pt) of event-by-event hy¬ 
dro. The qualitative features of both principal compo¬ 
nents are reproduced well by single-shot hydrodynamics, 
especially for the leading flow. It is particularly notable 
how the single-shot evolution reproduces the change of 
sign in Vy (pt)- However, in an important pt range, 
Pt ^ 1.2 GeV, the single-shot evolution misses the event- 
by-event curve for subleading flow by ~ 30%. 


It is useful to examine the time development of the sub¬ 
leading flow in the single-shot hydrodynamics. In Fig. 11, 
we present three snapshots of the subleading flow evo¬ 
lution. The color contours show the radial momentum 
density per rapidity, 

rT”” = r(e -I- p)u^u '', (39) 

as a function of proper time r. 

Shortly after the formation of the fireball, at r = 2.6 fm 
we observe negative triangular flow in Fig. 11(a). This 
flow is produced by the excess of material at small radii 
flowing into the “valleys” at larger radii [see Fig. 4(b)]. 
However, the radial flow has not developed yet, and 
therefore this phase of the evolution creates negative flow 
at small transverse momentum. After this stage, we see 
typical flow evolution of a triangular perturbation, i.e. 
the negative geometric eccentricity at small radii is trans¬ 
formed into positive triangular flow at large transverse 
momentum [see Figs. 11(b) and (c)]. The inner eccen¬ 
tricity dominates over the outer eccentricity at high pT 
because the radial flow has more time to develop before 
freeze-out, and because there is more material at small 
radii. 


VI. DISCUSSION 

This paper illustrates how principal component anal¬ 
ysis can be used to understand the physics encoded in 
the two particle correlation matrix of hydrodynamics. 
PGA is an economical way to summarize the factoriza¬ 
tion breaking in these correlations. More precisely, we 
found that the r matrix of correlation coefficients in hy¬ 
drodynamics, Eq. (1), is completely described by two 
principal components, V^^\pt) and V^'^\pt) as written 
in Eq. (2). Importantly, these components have a sim¬ 
ple physical interpretation—they are the hydrodynamic 
response to two statistically independent initial condi¬ 
tions in the fluctuating geometry. The leading principle 
component is the hydrodynamic response to the partic¬ 
ipant triangularity, while the subleading flow (which is 
uncorrelated with the leading flow) is the hydrodynamic 
response to the first radial excitation of the triangularity. 
This conclusion was reached by averaging the event-by- 
event geometry in the event plane of the subleading flow 
(Fig. 4). The magnitude of this radial excitation is on 
par with the magnitude of the triangularity (Fig. 5), al¬ 
though the hydro response is smaller in magnitude. Since 
the subleading component is uncorrelated with the inte¬ 
grated U 3 , it is projected out in analyses of triangular 
flow based on the scalar product or event plane methods. 

We first studied the basic properties of the sublead¬ 
ing flow such as its dependence on transverse momen¬ 
tum (Fig. 1), and centrality and shear viscosity (Fig. 3). 
The flow response is approximately linear to the geomet¬ 
rical deformation. This was checked by simulating the 
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r = 2.6fm r = 5.6fm r = 8.6fm 
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FIG. 11. Hydrodynamic evolution of the subleading triangular flow for the averaged initial conditions shown in Fig. 4(b). The 
color contours indicate the radial momentum density per rapidity, while the arrows indicate the radial flow velocity. 


response to the average in-plane geometry with “single¬ 
shot” hydrodynamics (Fig. 11), and comparing this re¬ 
sult to event-by-event hydrodynamics; i.e., we compared 
the response to the average with the averaged response 
(Fig. 10). 

Motivated by the linearity of the response, we con¬ 
structed a geometrical predictor for the subleading flow 
analogous to £ 3 , 3 . We first defined £ 3 (fc) as the fc—th 
Fourier mode of the event-by-event triangular geometry 
up to normalization.® Then we constructed a linear geo¬ 
metrical predictor for the leading and subleading flow an¬ 
gles and magnitudes based on two Fourier modes (Figs. 7 
and 12). Indeed, the subleading flow response is propor¬ 
tional to an event-by-event quantity which captures the 
radial excitation of the triangular geometry, 

J (fx s(to, x) e*®'^p(r), (40) 

where s{to, x) is the initial entropy distribution, and p(r) 
is an appropriate excited radial weight function. The two 
term Fourier fit to p(r) is tabulated in the appendix and 
graphed in Fig. 8 . The subleading flow probes the initial 
state geometry at higher wave numbers than the leading 
flow (Fig. 6 ). We found that the correlation between 
the flow and the Fourier components of the geometry is 
maximized for wave numbers away from zero, fci?rms 
4.0. Thus, the subleading flow provides a new test of 
viscous hydrodynamics and initial-state models. 

In peripheral collisions the correlation between the 
linear geometrical predictor and the flow is smaller. 
This suggests that nonlinear dynamics at large impact 
parameters couples the average elliptic geometry to the 


® The normalization is chosen so that limj,_>o esffc) = £3,3. 


harmonic perturbations [9, 12, 18]. The statistical tools 
such as PCA and related methods developed in this 
work can be used to clarify this complex hydrodynamic 
response. 
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Appendix A: Two term predictor 


Here we present the best fit results for the two wave 
number predictor, see eqs. (34), (36) and (37), 


— 43,1; 43,2 TZq — 3 Ri-nis ■ (-^^) 


Table II records the ratios of fit coefficients UJ 2 /UJ 1 for the 
leading and subleading predictors. 

In Fig. 12 we show the correlations between the flow 
and its predictor for both the angles and magnitudes. 
The subleading flow direction correlates well with the 
predictor, and there is reasonable correlation for the mag¬ 
nitude as well. 
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FIG. 12. Angle and magnitude correlations between the flow and and the two term predictor in central collisions. 


centrality 

leading 

0J2I0J1 

subleading 

0J2I0J1 

0-5% 

- 0.93 

-2.61 

5-10% 

-0.87 

-2.89 

10-15% 

-0.83 

-3.03 

15-20% 

-0.78 

-3.08 

20-25% 

-0.73 

-3.16 

25-30% 

-0.70 

-3.18 

30-35% 

-0.63 

-3.22 

35-40% 

-0.58 

-3.17 

40-45% 

-0.53 

-3.15 

45-50% 

-0.44 

-3.12 

50-55% 

-0.35 

-3.07 

55-60% 

-0.23 

-3.02 

60-65% 

-0.06 

-2.99 

65-70% 

0.08 

-2.88 


TABLE II. Table of two term predictor coefficients, eqs. (34), 
(36) and (37). 
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